#ifndef _MESH_
#define _MESH_

#include "Dofs.h"
#include "Domain.h"
#include "Grid.h"
#include "math.h"
#include <fstream>
#include <functional>
#include <map>
#include <set>
#include <sstream>
#include <tbb/task.h>
using namespace std;
/**
 * @brief
 * 定义域划分成网格Mesh,Mesh包含小网格Grid,每个Grid的类型由TemplateElement确定
 *
 * @tparam int DIM 维度
 */

vector<int> GridId1toId2(int i, int _segx);
int DofsId2toId1(int i, int j, int _segx, int _order);

template <unsigned int DIM> class Mesh {
  public:
    /**
     * @brief 默认构造函数
     */
    Mesh() = default;

    /**
     * @brief 默认析构函数
     */
    virtual ~Mesh() = default;
    /**
     * @brief 返回自由度的总个数
     * @details _GridList.size()*getNumDofs();
     * @return
     */
    virtual int getTotalNumDofs() = 0;
    /**
     * @brief 返回网格的总个数
     * @details _GridList.size()
     * @return [description]
     */
    virtual int getNumGrid() = 0;
    /**
     * @brief 返回存储Grid的vector
     * @details _GridList
     * @return [description]
     */
    virtual vector<Grid<DIM>> &getGrid() = 0;
    /**
     * @brief 返回序号为_index的Grid
     *
     * @param _index Grid的序号
     * @return _GridList[i];
     */
    virtual Grid<DIM> &getGrid(int _index) = 0;
    /**
     * @brief   返回vector，存储边界上的自由度
     * @details
     * 对所有grid,依据seg找出边界的grid,再确定对应边界的自由度下标是哪几个
     *
     * @param component 试图对边界分段，并标号
     * @return 存储边界上自由度的vector
     */
    virtual vector<Grid<2>> getBoundaryGrid() = 0;
    // virtual void setBoundaryIndex() = 0;
    virtual set<Dofs<2>> getBoundaryDofs() = 0;
    virtual vector<pair<int, int>> getIndexofGrid(int i) = 0;
    virtual vector<int> &getSegment() = 0;
};

class Q1Mesh : public Mesh<2> {
  public:
    Q1Mesh() = default;
    /**
     * @brief 含参构造函数,根据给定的定义域和划分段数来划分网格。
     * @param domain 定义域
     * @param seg 每个维度要划分的段数
     */
    Q1Mesh(Domain<2> *_domain, vector<int> _seg);
    ~Q1Mesh() = default;

    int getTotalNumDofs();
    int getNumGrid();
    vector<Grid<2>> &getGrid();
    Grid<2> &getGrid(int i);
    vector<Grid<2>> getBoundaryGrid();
    vector<int> &getSegment();
    // void setBoundaryIndex();
    set<Dofs<2>> getBoundaryDofs();
    vector<pair<int, int>> getIndexofGrid(int i);

  protected:
    vector<Grid<2>> gridList; /*< 存放所有的Grid*/
    vector<pair<int, int>>
        indexonMesh; /*< 存放mesh上的index,dofsindex,boundarymark */
    vector<int> seg; /*< 存放各维度分段数*/
};

Q1Mesh::Q1Mesh(Domain<2> *_domain, vector<int> _seg) {
    assert(_seg.size() == 2);
    Point<2> p0 = _domain->getVertex(0);
    Point<2> p2 = _domain->getVertex(2);
    double xLen = (p2[0] - p0[0]) / (1.0 * _seg[0]);
    double yLen = (p2[1] - p0[1]) / (1.0 * _seg[1]);
    for (int j = 0; j < _seg[1]; ++j) {
        for (int i = 0; i < _seg[0]; ++i) {
            vector<Point<2>> pnts;
            pnts.push_back({p0[0] + i * xLen, p0[1] + j * yLen});
            pnts.push_back({p0[0] + (i + 1) * xLen, p0[1] + j * yLen});
            pnts.push_back({p0[0] + (i + 1) * xLen, p0[1] + (j + 1) * yLen});
            pnts.push_back({p0[0] + i * xLen, p0[1] + (j + 1) * yLen});
            Grid<2> grid(pnts);
            gridList.push_back(grid);
        }
    }
    for (int i = 0; i < (_seg[0] + 1) * (_seg[1] + 1); i++) {
        indexonMesh.push_back({i, 0});
    }
    seg = _seg;
}

int Q1Mesh::getNumGrid() { return gridList.size(); }

Grid<2> &Q1Mesh::getGrid(int i) { return gridList[i]; }

vector<Grid<2>> &Q1Mesh::getGrid() { return gridList; }

vector<Grid<2>> Q1Mesh::getBoundaryGrid() {
    vector<Grid<2>> boundary_grid;
    int k = 0;
    for (int i = 0; i < seg[0] - 1; i++) {
        boundary_grid.push_back(gridList[k]);
        k++;
    }
    k = seg[0];
    for (int i = 0; i < seg[1] - 1; i++) {
        boundary_grid.push_back(gridList[k]);
        k += seg[0];
    }
    k = seg[0] - 1;
    for (int i = 0; i < seg[1] - 1; i++) {
        boundary_grid.push_back(gridList[k]);
        k += seg[0];
    }
    k = seg[0] * (seg[1] - 1);
    for (int i = 0; i < seg[0] - 1; i++) {
        boundary_grid.push_back(gridList[k]);
        k++;
    }
    return boundary_grid;
}

vector<int> &Q1Mesh::getSegment() { return seg; }

vector<pair<int, int>> Q1Mesh::getIndexofGrid(int i) {
    vector<pair<int, int>> idx;
    int a = seg[0];
    idx.push_back(indexonMesh[i + i / a]);
    idx.push_back(indexonMesh[i + i / a + 1]);
    idx.push_back(indexonMesh[i + i / a + a + 2]);
    idx.push_back(indexonMesh[i + i / a + a + 1]);
    return idx;
}

// void Q1Mesh::setBoundaryIndex()
// {
//     vector<int> k = {0, seg[0] + 1, seg[0], (seg[0] + 1) * seg[1] + 1};
//     for (int i = 0; i < seg[0]; i++)
//     {
//         indexonMesh[k[0]].second = 1;
//         k[0]++;
//         indexonMesh[k[3]].second = 1;
//         k[3]++;
//     }
//     for (int i = 0; i < seg[1]; i++)
//     {
//         indexonMesh[k[1]].second = 1;
//         k[1] += seg[0] + 1;
//         indexonMesh[k[2]].second = 1;
//         k[2] += seg[0] + 1;
//     }
// }

set<Dofs<2>> Q1Mesh::getBoundaryDofs() {
    set<Dofs<2>> boundaryDofs;
    vector<int> k = {0, 0, seg[0] - 1, seg[0] * (seg[1] - 1)};
    vector<int> m = {0, seg[0] + 1, seg[0], (seg[0] + 1) * seg[1] + 1};
    for (int i = 0; i < seg[0]; i++) {
        Dofs<2> d0 = {gridList[k[0]].getVertices(0), m[0]};
        Dofs<2> d3 = {gridList[k[3]].getVertices(2), m[3]};
        boundaryDofs.insert(d0);
        boundaryDofs.insert(d3);
        k[0]++;
        k[3]++;
        m[0]++;
        m[3]++;
    }
    for (int i = 0; i < seg[1]; i++) {
        Dofs<2> d1 = {gridList[k[1]].getVertices(3), m[1]};
        Dofs<2> d2 = {gridList[k[2]].getVertices(1), m[2]};
        boundaryDofs.insert(d1);
        boundaryDofs.insert(d2);
        k[1] += seg[0];
        k[2] += seg[0];
        m[1] += seg[0] + 1;
        m[2] += seg[0] + 1;
    }
    return boundaryDofs;
}

int Q1Mesh::getTotalNumDofs() { return (seg[0] + 1) * (seg[1] + 1); }

class Q2Mesh : public Mesh<2> {
  public:
    Q2Mesh() = default;
    /**
     * @brief 含参构造函数,根据给定的定义域和划分段数来划分网格。
     * @param domain 定义域
     * @param seg 每个维度要划分的段数
     */
    Q2Mesh(Domain<2> *_domain, vector<int> _seg);
    ~Q2Mesh() = default;

    int getTotalNumDofs();
    int getNumGrid();
    vector<Grid<2>> &getGrid();
    Grid<2> &getGrid(int i);
    vector<Grid<2>> getBoundaryGrid();
    vector<int> &getSegment();
    // void setBoundaryIndex();
    set<Dofs<2>> getBoundaryDofs();
    vector<pair<int, int>> getIndexofGrid(int i);

  protected:
    vector<Grid<2>> gridList;           /*< 存放所有的Grid*/
    vector<pair<int, int>> indexonMesh; /*< 存放mesh上的index */
    vector<int> seg;                    /*< 存放各维度分段数*/
};

Q2Mesh::Q2Mesh(Domain<2> *_domain, vector<int> _seg) {
    assert(_seg.size() == 2);
    Point<2> p0 = _domain->getVertex(0);
    Point<2> p2 = _domain->getVertex(2);
    double xLen = (p2[0] - p0[0]) / (1.0 * _seg[0]);
    double yLen = (p2[1] - p0[1]) / (1.0 * _seg[1]);
    for (int j = 0; j < _seg[1]; ++j) {
        for (int i = 0; i < _seg[0]; ++i) {
            vector<Point<2>> pnts;
            pnts.push_back({p0[0] + i * xLen, p0[1] + j * yLen});
            pnts.push_back({p0[0] + (i + 1) * xLen, p0[1] + j * yLen});
            pnts.push_back({p0[0] + (i + 1) * xLen, p0[1] + (j + 1) * yLen});
            pnts.push_back({p0[0] + i * xLen, p0[1] + (j + 1) * yLen});
            Grid<2> grid(pnts);
            gridList.push_back(grid);
        }
    }
    for (int i = 0; i < (2 * _seg[0] + 1) * (2 * _seg[1] + 1); i++) {
        indexonMesh.push_back({i, 0});
    }
    seg = _seg;
}

int Q2Mesh::getNumGrid() { return gridList.size(); }

Grid<2> &Q2Mesh::getGrid(int i) { return gridList[i]; }

vector<Grid<2>> &Q2Mesh::getGrid() { return gridList; }

vector<int> &Q2Mesh::getSegment() { return seg; }

int Q2Mesh::getTotalNumDofs() { return (2 * seg[0] + 1) * (2 * seg[1] + 1); }

// void Q2Mesh::setBoundaryIndex()
// {
//     vector<int> k = {0, 2 * seg[0] + 1, 2 * seg[0], (2 * seg[0] + 1) * (2 *
//     seg[1]) + 1}; for (int i = 0; i < 2 * seg[0]; i++)
//     {
//         indexonMesh[k[0]].second = 1;
//         k[0]++;
//         indexonMesh[k[3]].second = 1;
//         k[3]++;
//     }
//     for (int i = 0; i < 2 * seg[1]; i++)
//     {
//         indexonMesh[k[1]].second = 1;
//         k[1] += 2 * seg[0] + 1;
//         indexonMesh[k[2]].second = 1;
//         k[2] += 2 * seg[0] + 1;
//     }
// }

vector<pair<int, int>> Q2Mesh::getIndexofGrid(int i) {
    vector<pair<int, int>> idx;
    int a = seg[0];
    idx.push_back(indexonMesh[i / a * (4 * a + 2) + i % a * 2]);
    idx.push_back(indexonMesh[i / a * (4 * a + 2) + i % a * 2 + 1]);
    idx.push_back(indexonMesh[i / a * (4 * a + 2) + i % a * 2 + 2]);
    idx.push_back(indexonMesh[i / a * (4 * a + 2) + 2 * a + 1 + i % a * 2 + 2]);
    idx.push_back(indexonMesh[(i / a + 1) * (4 * a + 2) + i % a * 2 + 2]);
    idx.push_back(indexonMesh[(i / a + 1) * (4 * a + 2) + i % a * 2 + 1]);
    idx.push_back(indexonMesh[(i / a + 1) * (4 * a + 2) + i % a * 2]);
    idx.push_back(indexonMesh[i / a * (4 * a + 2) + 2 * a + 1 + i % a * 2]);
    idx.push_back(indexonMesh[i / a * (4 * a + 2) + 2 * a + 1 + i % a * 2 + 1]);
    return idx;
}

vector<Grid<2>> Q2Mesh::getBoundaryGrid() {
    vector<Grid<2>> boundary_grid;
    int k = 0;
    for (int i = 0; i < seg[0] - 1; i++) {
        boundary_grid.push_back(gridList[k]);
        k++;
    }
    k = seg[0];
    for (int i = 0; i < seg[1] - 1; i++) {
        boundary_grid.push_back(gridList[k]);
        k += seg[0];
    }
    k = seg[0] - 1;
    for (int i = 0; i < seg[1] - 1; i++) {
        boundary_grid.push_back(gridList[k]);
        k += seg[0];
    }
    k = seg[0] * (seg[1] - 1);
    for (int i = 0; i < seg[0] - 1; i++) {
        boundary_grid.push_back(gridList[k]);
        k++;
    }
    return boundary_grid;
}

set<Dofs<2>> Q2Mesh::getBoundaryDofs() {
    set<Dofs<2>> boundaryDofs;
    vector<int> k = {0, 0, seg[0] - 1, seg[0] * (seg[1] - 1)};
    vector<int> m = {0, 2 * seg[0] + 1, 2 * seg[0],
                     (2 * seg[0] + 1) * 2 * seg[1] + 1};
    for (int i = 0; i < seg[0]; i++) {
        Point<2> p0_0 = gridList[k[0]].getVertices(0);
        Point<2> p0_1 = gridList[k[0]].getVertices(1);
        Dofs<2> d0_0 = {p0_0, m[0]};
        Dofs<2> d0_1 = {(p0_0 + p0_1) / 2, m[0] + 1};
        Point<2> p2_2 = gridList[k[3]].getVertices(2);
        Point<2> p2_3 = gridList[k[3]].getVertices(3);
        Dofs<2> d2_5 = {(p2_2 + p2_3) / 2, m[3]};
        Dofs<2> d2_4 = {p2_2, m[3] + 1};
        boundaryDofs.insert(d0_0);
        boundaryDofs.insert(d0_1);
        boundaryDofs.insert(d2_5);
        boundaryDofs.insert(d2_4);
        k[0]++;
        k[3]++;
        m[0] += 2;
        m[3] += 2;
    }
    for (int i = 0; i < seg[1]; i++) {
        Point<2> p3_0 = gridList[k[1]].getVertices(0);
        Point<2> p3_3 = gridList[k[1]].getVertices(3);
        Dofs<2> d3_7 = {(p3_0 + p3_3) / 2, m[1]};
        Dofs<2> d3_6 = {p3_3, m[1] + 2 * seg[0] + 1};
        Point<2> p1_1 = gridList[k[2]].getVertices(1);
        Point<2> p1_2 = gridList[k[2]].getVertices(2);
        Dofs<2> d1_2 = {p1_1, m[2]};
        Dofs<2> d1_3 = {(p1_1 + p1_2) / 2, m[2] + 2 * seg[0] + 1};
        boundaryDofs.insert(d3_7);
        boundaryDofs.insert(d3_6);
        boundaryDofs.insert(d1_2);
        boundaryDofs.insert(d1_3);
        k[1] += seg[0];
        k[2] += seg[0];
        m[1] += 2 * (2 * seg[0] + 1);
        m[2] += 2 * (2 * seg[0] + 1);
    }
    return boundaryDofs;
}

class Q4Mesh : public Mesh<2> {
  public:
    Q4Mesh() = default;
    /**
     * @brief 含参构造函数,根据给定的定义域和划分段数来划分网格。
     * @param domain 定义域
     * @param seg 每个维度要划分的段数
     */
    Q4Mesh(Domain<2> *_domain, vector<int> _seg);
    ~Q4Mesh() = default;

    int getTotalNumDofs();
    int getNumGrid();
    vector<Grid<2>> &getGrid();
    Grid<2> &getGrid(int i);
    vector<Grid<2>> getBoundaryGrid();
    vector<int> &getSegment();
    // void setBoundaryIndex();
    set<Dofs<2>> getBoundaryDofs();
    vector<pair<int, int>> getIndexofGrid(int i);

  protected:
    vector<Grid<2>> gridList;           /*< 存放所有的Grid*/
    vector<pair<int, int>> indexonMesh; /*< 存放mesh上的index */
    vector<int> seg;                    /*< 存放各维度分段数*/
};

Q4Mesh::Q4Mesh(Domain<2> *_domain, vector<int> _seg) {
    assert(_seg.size() == 2);
    Point<2> p0 = _domain->getVertex(0);
    Point<2> p2 = _domain->getVertex(2);
    double xLen = (p2[0] - p0[0]) / (1.0 * _seg[0]);
    double yLen = (p2[1] - p0[1]) / (1.0 * _seg[1]);
    for (int j = 0; j < _seg[1]; ++j) {
        for (int i = 0; i < _seg[0]; ++i) {
            vector<Point<2>> pnts;
            pnts.push_back({p0[0] + i * xLen, p0[1] + j * yLen});
            pnts.push_back({p0[0] + (i + 1) * xLen, p0[1] + j * yLen});
            pnts.push_back({p0[0] + (i + 1) * xLen, p0[1] + (j + 1) * yLen});
            pnts.push_back({p0[0] + i * xLen, p0[1] + (j + 1) * yLen});
            Grid<2> grid(pnts);
            gridList.push_back(grid);
        }
    }
    for (int i = 0; i < (4 * _seg[0] + 1) * (4 * _seg[1] + 1); i++) {
        indexonMesh.push_back({i, 0});
    }
    seg = _seg;
}

int Q4Mesh::getNumGrid() { return gridList.size(); }

Grid<2> &Q4Mesh::getGrid(int i) { return gridList[i]; }

vector<Grid<2>> &Q4Mesh::getGrid() { return gridList; }

vector<int> &Q4Mesh::getSegment() { return seg; }

int Q4Mesh::getTotalNumDofs() { return (4 * seg[0] + 1) * (4 * seg[1] + 1); }

vector<pair<int, int>> Q4Mesh::getIndexofGrid(int i) {
    vector<pair<int, int>> idx;
    int x = GridId1toId2(i, seg[0])[0];
    int y = GridId1toId2(i, seg[0])[1];
    for (int j = 0; j < 5; j++)
        for (int i = 0; i < 5; i++)
            idx.push_back(
                indexonMesh[DofsId2toId1(4 * x + i, 4 * y + j, seg[0], 4)]);
    return idx;
}

vector<Grid<2>> Q4Mesh::getBoundaryGrid() {
    vector<Grid<2>> boundary_grid;
    int k = 0;
    for (int i = 0; i < seg[0] - 1; i++) {
        boundary_grid.push_back(gridList[k]);
        k++;
    }
    k = seg[0];
    for (int i = 0; i < seg[1] - 1; i++) {
        boundary_grid.push_back(gridList[k]);
        k += seg[0];
    }
    k = seg[0] - 1;
    for (int i = 0; i < seg[1] - 1; i++) {
        boundary_grid.push_back(gridList[k]);
        k += seg[0];
    }
    k = seg[0] * (seg[1] - 1);
    for (int i = 0; i < seg[0] - 1; i++) {
        boundary_grid.push_back(gridList[k]);
        k++;
    }
    return boundary_grid;
}

set<Dofs<2>> Q4Mesh::getBoundaryDofs() {
    set<Dofs<2>> boundaryDofs;
    vector<int> k = {0, 0, seg[0] - 1, seg[0] * (seg[1] - 1)};
    vector<int> m = {0, 4 * seg[0] + 1, 4 * seg[0],
                     (4 * seg[0] + 1) * 4 * seg[1] + 1};
    for (int i = 0; i < seg[0]; i++) {
        Point<2> p0_0 = gridList[k[0]].getVertices(0);
        Point<2> p0_1 = gridList[k[0]].getVertices(1);
        Dofs<2> d0_0 = {p0_0, m[0]};
        Dofs<2> d0_1 = {p0_0 * 0.75 + p0_1 * 0.25, m[0] + 1};
        Dofs<2> d0_2 = {p0_0 * 0.5 + p0_1 * 0.5, m[0] + 2};
        Dofs<2> d0_3 = {p0_0 * 0.25 + p0_1 * 0.75, m[0] + 3};
        Point<2> p2_2 = gridList[k[3]].getVertices(2);
        Point<2> p2_3 = gridList[k[3]].getVertices(3);
        Dofs<2> d2_21 = {p2_3 * 0.75 + p2_2 * 0.25, m[3]};
        Dofs<2> d2_22 = {p2_3 * 0.5 + p2_2 * 0.5, m[3] + 1};
        Dofs<2> d2_23 = {p2_3 * 0.25 + p2_2 * 0.75, m[3] + 2};
        Dofs<2> d2_24 = {p2_2, m[3] + 3};
        boundaryDofs.insert(d0_0);
        boundaryDofs.insert(d0_1);
        boundaryDofs.insert(d0_2);
        boundaryDofs.insert(d0_3);
        boundaryDofs.insert(d2_21);
        boundaryDofs.insert(d2_22);
        boundaryDofs.insert(d2_23);
        boundaryDofs.insert(d2_24);
        k[0]++;
        k[3]++;
        m[0] += 4;
        m[3] += 4;
    }
    for (int i = 0; i < seg[1]; i++) {
        Point<2> p3_0 = gridList[k[1]].getVertices(0);
        Point<2> p3_3 = gridList[k[1]].getVertices(3);
        Dofs<2> d3_5 = {p3_0 * 0.75 + p3_3 * 0.25, m[1]};
        Dofs<2> d3_10 = {p3_0 * 0.5 + p3_3 * 0.5, m[1] + 4 * seg[0] + 1};
        Dofs<2> d3_15 = {p3_0 * 0.25 + p3_3 * 0.75,
                         m[1] + 2 * (4 * seg[0] + 1)};
        Dofs<2> d3_20 = {p3_3, m[1] + 3 * (4 * seg[0] + 1)};
        Point<2> p1_1 = gridList[k[2]].getVertices(1);
        Point<2> p1_2 = gridList[k[2]].getVertices(2);
        Dofs<2> d1_4 = {p1_1, m[2]};
        Dofs<2> d1_9 = {p1_1 * 0.75 + p1_2 * 0.25, m[2] + 4 * seg[0] + 1};
        Dofs<2> d1_14 = {p1_1 * 0.5 + p1_2 * 0.5, m[2] + 2 * (4 * seg[0] + 1)};
        Dofs<2> d1_19 = {p1_1 * 0.25 + p1_2 * 0.75,
                         m[2] + 3 * (4 * seg[0] + 1)};
        boundaryDofs.insert(d3_5);
        boundaryDofs.insert(d3_10);
        boundaryDofs.insert(d3_15);
        boundaryDofs.insert(d3_20);
        boundaryDofs.insert(d1_4);
        boundaryDofs.insert(d1_9);
        boundaryDofs.insert(d1_14);
        boundaryDofs.insert(d1_19);
        k[1] += seg[0];
        k[2] += seg[0];
        m[1] += 4 * (4 * seg[0] + 1);
        m[2] += 4 * (4 * seg[0] + 1);
    }
    return boundaryDofs;
}

vector<int> GridId1toId2(int i, int _segx) {
    vector<int> a = {i % _segx, i / _segx};
    return a;
}

int DofsId2toId1(int i, int j, int _segx, int _order) {
    return i + (_order * _segx + 1) * j;
}



// class EasyMesh : public Mesh<2> {
//   private:
//     string meshName;
//     long int nDofs;
//     long int nGrids;
//     long int nEdges;
//     long int nBoundaryEdges;
//     vector<Dofs<2>> dofsList;
//     vector<Dofs<2>> boundaryDofsList;
//     vector<Grid<2>> gridsList;       /*< 存放所有的Grid*/
//     vector<vector<int>> indexonGrid; /*存放每个grid上的dofsindex*/
//     vector<pair<int, int>>
//         indexonMesh; /*存放所有的自由度编号和boundarymark对,其实编号和下标是一样的……*/
//     vector<vector<int>> boundaryEdgesList;
//     // not use function
//     vector<int> seg = {-1, -1};
//     virtual vector<int> &getSegment() override { return seg; };

//   public:
//     EasyMesh() = default;
//     EasyMesh(string meshname);
//     void inputDofData();
//     void inputEdgeData();
//     void inputGridData();
//     virtual int getNumGrid() override { return nGrids; };
//     virtual int getTotalNumDofs() override { return nDofs; };
//     virtual Grid<2> &getGrid(int index) override;
//     virtual vector<Grid<2>> &getGrid() override;
//     virtual vector<Grid<2>> getBoundaryGrid() override;
//     virtual set<Dofs<2>> getBoundaryDofs() override;
//     virtual vector<pair<int, int>> getIndexofGrid(int i) override;
//     int getBoundaryMark(int index);
//     vector<Dofs<2>> &getAllDofs();
// };

// vector<Dofs<2>> &EasyMesh::getAllDofs() { return dofsList; }
// int EasyMesh::getBoundaryMark(int index) {
//     pair<int, int> index_boundarymark = indexonMesh.at(index);
//     return index_boundarymark.second;
// }
// EasyMesh::EasyMesh(string meshname) {
//     meshName = meshname;
//     inputDofData();
//     inputEdgeData();
//     inputGridData();
// }

// void EasyMesh::inputDofData() {
//     ifstream fin;
//     fin.open(meshName + ".n");
//     string line;
//     stringstream ss;
//     getline(fin, line);
//     ss.clear();
//     ss.str(line);
//     ss >> nDofs >> nGrids >> nEdges;
//     int index, bndflag;
//     double x, y;
//     for (int i = 0; i < nDofs; i++) {
//         getline(fin, line);
//         ss.clear();
//         ss.str(line);
//         ss >> index >> x >> y >> bndflag;

//         Point<2> temp_point({x, y});
//         Dofs<2> temp_dof(temp_point, index);
//         indexonMesh.push_back({index, bndflag});
//         dofsList.push_back(temp_dof);
//         if (bndflag)
//             boundaryDofsList.push_back(temp_dof);
//     }
// }

// void EasyMesh::inputEdgeData() {
//     ifstream fin;
//     fin.open(meshName + ".s");
//     string line;
//     stringstream ss;
//     getline(fin, line);
//     ss.clear();
//     ss.str(line);
//     // int _nedge;
//     // ss >> _nedge;
//     int a, d0, d1, bndflag;
//     for (int i = 0; i < nEdges; i++) {
//         getline(fin, line);
//         ss.clear();
//         ss.str(line);
//         ss >> a >> d0 >> d1 >> a >> a >> bndflag;
//         if (bndflag) {
//             vector<int> temp_edge(2);
//             temp_edge[0] = d0;
//             temp_edge[1] = d1;
//             boundaryEdgesList.push_back(temp_edge);
//         }
//     }
//     nBoundaryEdges = boundaryEdgesList.size();
// }

// void EasyMesh::inputGridData() {
//     ifstream fin;
//     fin.open(meshName + ".e");
//     string line;
//     stringstream ss;
//     getline(fin, line);

//     int a, d0, d1, d2;
//     for (int i = 0; i < nGrids; i++) {
//         getline(fin, line);
//         ss.clear();
//         ss.str(line);
//         ss >> a >> d0 >> d1 >> d2;
//         vector<int> indexs = {d0, d1, d2};
//         indexonGrid.push_back(indexs);
//         vector<Point<2>> points = {dofsList[d0].getPoint(),
//                                    dofsList[d1].getPoint(),
//                                    dofsList[d2].getPoint()};
//         Grid<2> grid(points);
//         gridsList.push_back(grid);
//     }
// };

// Grid<2> &EasyMesh::getGrid(int index) { return gridsList.at(index); };
// vector<Grid<2>> &EasyMesh::getGrid() { return gridsList; };
// vector<Grid<2>> EasyMesh::getBoundaryGrid() {
//     // why we need to getBoundaryGrid? we already have getBoundaryDofs
//     // TODO:this function is not used,we can delete it.
//     return gridsList;
// };

// set<Dofs<2>> EasyMesh::getBoundaryDofs() {
//     // use set to avoid duplicate dofs
//     // unnecessary in unstructed mesh
//     set<Dofs<2>> dofsSet(boundaryDofsList.begin(), boundaryDofsList.end());
//     return dofsSet;
// };

// vector<pair<int, int>> EasyMesh::getIndexofGrid(int i) {
//     vector<pair<int, int>> res;
//     vector<int> dofsindex = indexonGrid.at(i);
//     for (int &index : dofsindex) {
//         auto boundarymark = find_if(
//             indexonMesh.begin(), indexonMesh.end(),
//             [index](const pair<int, int> &p) { return p.first == index; });
//         if (boundarymark != indexonMesh.end()) {
//             res.emplace_back(index, (*boundarymark).second);
//         } else {
//             std::cerr << "no such index in indexonMesh" << std::endl;
//         }
//     }
//     return res;
// };

// TODO:从easymesh中生成三次元网格。
class EasyMesh : public Mesh<2> {
    using elementIdx = int;
    using edgeIdx = int;

  private:
    unsigned int Order = 1;
    string meshName;
    long int nDofs;                   //需要重新分配
    long int nGrids;                  //不变
    long int nEdges;                  //不变
    long int nBoundaryEdges;          //不变
    vector<Dofs<2>> dofsList;         /*<存放所有的自由度*/
    vector<Dofs<2>> boundaryDofsList; /*<存放所有边界上的自由度*/
    vector<Grid<2>> gridsList;        /*< 存放所有的Grid*/
    vector<vector<int>> indexonGrid;  /*存放每个grid上的dofsindex*/
    vector<vector<pair<elementIdx, edgeIdx>>>
        gridHasEdge;                    /*<edge是哪几个grid的几号边*/
    vector<vector<edgeIdx>> edgeofGrid; /*一个grid的几条边的编号*/
    vector<pair<int, int>>
        indexonMesh; /*存放所有的自由度编号和boundarymark对,其实编号和下标是一样的……*/
    vector<vector<int>> edgesList; /*<所有的边界的顶点的编号*/
    vector<vector<int>> boundaryEdgesList; /*< 存放所有边界上的edge的顶点的编号*/

    void inputDofData();
    void inputEdgeData();
    void inputGridData();
 
    /**
     * @brief
     * 根据当前grid的edge上的自由度的局部编号，推断出在该自由度在相邻grid中的局部编号
     *
     * @param[in]  edge1_index  edge在当前grid中的编号
     * @param[in]  local_index  dof在当前grid中的局部编号
     * @param[in]  edge2_index  edge在另一个grid中的编号
     *
     * @return     dof在另一个grid中的局部编号
     */
    int indexBetweenElement(int edge1_index, int local_index, int edge2_index);

    // not use function
    vector<int> seg = {-1, -1};
    virtual vector<int> &getSegment() override { return seg; };

  public:
    EasyMesh() = default;
    EasyMesh(string meshname,unsigned int order);
    virtual int getNumGrid() override { return nGrids; };
    virtual int getTotalNumDofs() override { return nDofs; };
    virtual Grid<2> &getGrid(int index) override;
    virtual vector<Grid<2>> &getGrid() override;
    virtual vector<Grid<2>> getBoundaryGrid() override;
    virtual set<Dofs<2>> getBoundaryDofs() override;
    virtual vector<pair<int, int>> getIndexofGrid(int i) override;
    int getBoundaryMark(int index);
    vector<Dofs<2>> &getAllDofs();
    vector<vector<Dofs<2> > > getAllEdges();
    int getNumEdge(){return nEdges;};
    Dofs<2>& getDofs(int index);
    void refine();
    void distributeDofs();

};

Dofs<2>& EasyMesh::getDofs(int index)
{
    assert(index >=0 && index < nDofs);
    return dofsList.at(index);
};

vector<Dofs<2>> &EasyMesh::getAllDofs() { return dofsList; }
int EasyMesh::getBoundaryMark(int index) {
    pair<int, int> index_boundarymark = indexonMesh.at(index);
    return index_boundarymark.second;
}
vector<vector<Dofs<2> > > EasyMesh::getAllEdges()
{
  // cout << "正在获取所有边:\n";
    vector<vector<Dofs<2>>> res;
    for (int edge_idx = 0; edge_idx < nEdges; edge_idx++){
        vector<int> indexs  =  edgesList.at(edge_idx);
        // cout << "edge_idx: " << edge_idx << "\n";
        // cout << "size of indexs:" << indexs.size() <<"\n";
        int start = indexs.at(0);
        int end = indexs.at(1);
        // cout << "start: " << start << ", end: " << end << "\n";
        vector<Dofs<2>> dofs_edge = {dofsList.at(start),dofsList.at(end)};
        res.push_back(dofs_edge);
    }
    // cout << "size of res edges: " << res.size() <<endl;
    return res;
};



EasyMesh::EasyMesh(string meshname,unsigned int order) {
    Order = order;
    meshName = meshname;
    inputDofData();
    inputEdgeData();
    inputGridData();
    // distributeDofs();
}

int EasyMesh::indexBetweenElement(int edge1_index, int local_index,
                                    int edge2_index) {
    vector<int> edge0_sec;
    for(int i = 0; i < Order - 1;++i) {
      edge0_sec.push_back(i+3);
    }
    vector<int> edge1_sec;
    for (auto dof : edge0_sec) {
        edge1_sec.push_back(dof + Order - 1);
    }
    vector<int> edge2_sec;
    for (auto dof : edge1_sec) {
        edge2_sec.insert(edge2_sec.begin(), dof + Order - 1);
    }
    vector<vector<int>> edge_sec = {edge0_sec, edge1_sec, edge2_sec};
    auto sec1 = edge_sec.at(edge1_index);
    auto sec2 = edge_sec.at(edge2_index);
    for (int i = 0; i < sec1.size(); ++i) {
        if (local_index == sec1.at(i)) {
            return sec2.at(i);
        }
    }
    std::cerr << "indexBetweenElement fails to get local index.\n";
    exit(EXIT_FAILURE);
};

void EasyMesh::refine() {
  cout << "refine the Mesh: " <<"\n";
  //重新构建所有的网格信息
  int new_nDofs = nDofs + nEdges*1;
  int new_nGrids = 4*nGrids;
  int new_nEdges = 2*nEdges + 3*nGrids;
  int new_nBoundaryEdges = 2* nBoundaryEdges;
   vector<Dofs<2>> new_dofsList =  dofsList;         /*<存放所有的自由度*/
    vector<Dofs<2>> new_boundaryDofsList = boundaryDofsList; /*<存放所有边界上的自由度*/
    vector<Grid<2>> new_gridsList(new_nGrids);        /*< 存放所有的Grid*/
    vector<vector<int>> new_indexonGrid(new_nGrids);// = indexonGrid;  /*存放每个grid上的dofsindex*/
    vector<vector<pair<elementIdx, edgeIdx>>>
        new_gridHasEdge(new_nEdges);                    /*<edge是哪几个grid的几号边*/
    vector<vector<edgeIdx>> new_edgeofGrid(new_nGrids); /*一个grid的几条边的编号*/
    vector<pair<int, int>>
        new_indexonMesh = indexonMesh; /*存放所有的自由度编号和boundarymark对,其实编号和下标是一样的……*/
    vector<vector<int>> new_edgesList(new_nEdges); /*<所有的边界的顶点的编号*/
    vector<vector<int>> new_boundaryEdgesList; /*< 存放所有边界上的edge的顶点的编号*/
    //中间信息:每个旧grid含有9条新边，6个新自由度
    vector<vector<int>> new_edges_old_grid = edgeofGrid;
    for(auto& edges_idx : new_edges_old_grid){
      edges_idx.resize(9);
    }
    vector<vector<int>> new_dofs_old_grid = indexonGrid;
    for(auto& dofs_idx : new_dofs_old_grid) {
      dofs_idx.resize(6);
    }
    int cglobal_dof_index = nDofs;
    for(int edge_idx = 0; edge_idx < nEdges; edge_idx++) {
      //取中点，构建新自由度，加入dofsList
      vector<int> end_pnt = edgesList.at(edge_idx);
      Dofs<2> start = dofsList.at(end_pnt.at(0));
      Dofs<2> end = dofsList.at(end_pnt.at(1));
      Point<2> new_pnt = (start.getPoint() + end.getPoint())/2.0;
      Dofs<2> new_dof(new_pnt,cglobal_dof_index);
      new_dofsList.push_back(new_dof);
      //新自由度的边界标识符确定，判断是否加入boundaryDofsList
      int start_bndmark = getBoundaryMark(end_pnt.at(0));
      int end_bndmark = getBoundaryMark(end_pnt.at(1));
      int new_bndmark = (start_bndmark==0||end_bndmark==0)?0:start_bndmark;
      if(new_bndmark != 0) {
        new_boundaryDofsList.push_back(new_dof);
      }
      //将index-boundarymark组成的pair加入new_indexonMesh
      pair<int,int> idx_bndmark = make_pair(cglobal_dof_index,new_bndmark);
      new_indexonMesh.push_back(idx_bndmark);
      //构建出两条新的边,加入new_edgesList
      vector<int> new_edge1 = {end_pnt.at(0),cglobal_dof_index};
      int new_edge_idx1 = 2*edge_idx;
      vector<int> new_edge2 = {cglobal_dof_index,end_pnt.at(1)};
      int new_edge_idx2 = 2*edge_idx+1;
        new_edgesList.at(new_edge_idx1) = new_edge1;
        new_edgesList.at(new_edge_idx2) = new_edge2;
      if(new_bndmark != 0) {
        if(start_bndmark != 0) {
          new_boundaryEdgesList.push_back(new_edge1);
        }
        if(end_bndmark != 0) {
          new_boundaryEdgesList.push_back(new_edge2);
        }
      }
      //根据gridHasEdge中这条边所在的单元，在new_edges_old_grid和new_dofs_old_grid写入对应位置。
      vector<pair<elementIdx, edgeIdx>> ee = gridHasEdge.at(edge_idx);
      for(auto& ee_idx : ee) {
        int gridIdx = ee_idx.first;
        int local_edgeIdx = ee_idx.second;
        //自由度对应0-3,1-4,2-5
        int local_dofIdx = local_edgeIdx + 3;
        new_dofs_old_grid.at(gridIdx).at(local_dofIdx) = cglobal_dof_index;
        int local_edge_idx1 = 2*local_edgeIdx;
        int local_edge_idx2 = 2*local_edgeIdx + 1;
        new_edges_old_grid.at(gridIdx).at(local_edge_idx1) = new_edge_idx1;
        new_edges_old_grid.at(gridIdx).at(local_edge_idx2) = new_edge_idx2;
      }

      cglobal_dof_index++;
    }
    //分配每个grid内部的三条新边，加入new_edgesList和new_edges_old_grid;
    int cglobal_edge_idx = 2 * nEdges;
    for(int grid_idx = 0; grid_idx < nGrids; grid_idx++) {
      vector<int> dofs_idx = new_dofs_old_grid.at(grid_idx);
      vector<int> edge6 = {dofs_idx.at(3),dofs_idx.at(5)};
      vector<int> edge7 = {dofs_idx.at(3),dofs_idx.at(4)};
      vector<int> edge8 = {dofs_idx.at(5),dofs_idx.at(4)};
      new_edgesList.at(cglobal_edge_idx) = edge6; 
      new_edges_old_grid.at(grid_idx).at(6) = cglobal_edge_idx;
      cglobal_edge_idx++;
      new_edgesList.at(cglobal_edge_idx) = edge7; 
      new_edges_old_grid.at(grid_idx).at(7) = cglobal_edge_idx;
      cglobal_edge_idx++;
      new_edgesList.at(cglobal_edge_idx) = edge8; 
      new_edges_old_grid.at(grid_idx).at(8) = cglobal_edge_idx;
      cglobal_edge_idx++;
    }
    // cout << "cglobal edge idx is:" << cglobal_edge_idx<< ", new_nEdges is :" << new_nEdges << "\n";
    //根据new_edges_old_grid和new_dofs_old_grid构建新grid,加入gridsList
    //构建indexonGrid和edgeofGrid
    for(int grid_idx = 0; grid_idx < nGrids; grid_idx++) {
      vector<int> dofs_idx = new_dofs_old_grid.at(grid_idx);
      vector<Dofs<2>> dofs;
      vector<Point<2>> pnts;
      for(auto idx : dofs_idx) {
        Dofs<2> dof = new_dofsList.at(idx);
        dofs.push_back(dof);
        pnts.push_back(dof.getPoint());
      }
      vector<int> edges_idx = new_edges_old_grid.at(grid_idx); 
      //0号,点:0-5-4,边:8-3-4
      vector<int> dofs_idx0 = {dofs_idx.at(0),dofs_idx.at(5),dofs_idx.at(4)};
      vector<Point<2>> pnts0 = {pnts.at(0),pnts.at(5),pnts.at(4)};
      vector<int> edges0 = {edges_idx.at(8),edges_idx.at(3),edges_idx.at(4)};
      Grid<2> grid0(pnts0);
      int grid_idx0 = 4*grid_idx;
      new_gridsList.at(grid_idx0) = grid0;
      new_indexonGrid.at(grid_idx0) = dofs_idx0;
      new_edgeofGrid.at(grid_idx0) = edges0;
      //1号,点:5-1-3，边:0-6-5
      vector<int> dofs_idx1 = {dofs_idx.at(5),dofs_idx.at(1),dofs_idx.at(3)};
      vector<Point<2>> pnts1 = {pnts.at(5),pnts.at(1),pnts.at(3)};
      vector<int> edges1 = {edges_idx.at(0),edges_idx.at(6),edges_idx.at(5)};
      Grid<2> grid1(pnts1);
      int grid_idx1 = 4*grid_idx+1;
      new_gridsList.at(grid_idx1) = grid1;
      new_indexonGrid.at(grid_idx1) = dofs_idx1;
      new_edgeofGrid.at(grid_idx1) = edges1;
      //2号,点:4-3-2,边:1-2-7
      vector<int> dofs_idx2 = {dofs_idx.at(4),dofs_idx.at(3),dofs_idx.at(2)};
      vector<Point<2>> pnts2 = {pnts.at(4),pnts.at(3),pnts.at(2)};
      vector<int> edges2 = {edges_idx.at(1),edges_idx.at(2),edges_idx.at(7)};
      Grid<2> grid2(pnts2);
      int grid_idx2 = 4*grid_idx+2;
      new_gridsList.at(grid_idx2) = grid2;
      new_indexonGrid.at(grid_idx2) = dofs_idx2;
      new_edgeofGrid.at(grid_idx2) = edges2;
      //3号,点:4-3-5，边:6-8-7
      vector<int> dofs_idx3 = {dofs_idx.at(4),dofs_idx.at(3),dofs_idx.at(5)};
      vector<Point<2>> pnts3 = {pnts.at(4),pnts.at(3),pnts.at(5)};
      vector<int> edges3 = {edges_idx.at(6),edges_idx.at(8),edges_idx.at(7)};
      Grid<2> grid3(pnts3);
      int grid_idx3 = 4*grid_idx+3;
      new_gridsList.at(grid_idx3) = grid3;
      new_indexonGrid.at(grid_idx3) = dofs_idx3;
      new_edgeofGrid.at(grid_idx3) = edges3;
    } 
    //根据edgeofGrid,构建gridHasEdge;
    for(int grid_idx = 0; grid_idx < new_nGrids; grid_idx++) {
      vector<int> edges_idx = new_edgeofGrid.at(grid_idx);
        auto p0 = make_pair(grid_idx, 0);
        new_gridHasEdge.at(edges_idx.at(0)).push_back(p0);
        auto p1 = make_pair(grid_idx, 1);
        new_gridHasEdge.at(edges_idx.at(1)).push_back(p1);
        auto p2 = make_pair(grid_idx, 2);
        new_gridHasEdge.at(edges_idx.at(2)).push_back(p2);
    }
    
    //修改所有数据
    nDofs = new_nDofs;
    nGrids = new_nGrids;
    nEdges = new_nEdges;
    nBoundaryEdges = new_nBoundaryEdges;
    dofsList = new_dofsList;         /*<存放所有的自由度*/
    boundaryDofsList = new_boundaryDofsList; /*<存放所有边界上的自由度*/
    gridsList =  new_gridsList;     /*< 存放所有的Grid*/
    indexonGrid = new_indexonGrid;  /*存放每个grid上的dofsindex*/
    gridHasEdge = new_gridHasEdge;  /*<edge是哪几个grid的几号边*/
    edgeofGrid = new_edgeofGrid; /*一个grid的几条边的编号*/
    indexonMesh = new_indexonMesh; /*存放所有的自由度编号和boundarymark*/
    edgesList = new_edgesList; /*<所有的边界的顶点的编号*/
    boundaryEdgesList = new_boundaryEdgesList; /*< 存放所有边界上的edge的顶点的编号*/
    cout << "自由度数量: " << nDofs << "\n";
    cout << "边数量: " << nEdges <<"\n";
    cout << "单元数量: " << nGrids << "\n";
    cout << "dofsList的大小: " << dofsList.size() <<"\n";
    cout << "gridsList的大小: " << gridsList.size() <<"\n";
    cout << "edgesList的大小: " << edgesList.size() <<"\n";

    // //如果依照grid循环，无法避免自由度的重复。
    // for(int gridIdx = 0; gridIdx < nGrids; gridIdx++)
    // {
    //   //构建新的自由度：边的中点，局部编号与二阶三角形单元相同。
    //   Grid<2> grid = gridsList.at(gridIdx);
    //   vector<int> dofs_idx = indexonGrid.at(gridIdx);
    //   vector<Point<2>> points = grid.getVertices();
    //   vector<int> bndmarks(3,0);
    //   for(int i = 0; i < 3; i++) {
    //     bndmarks.at(i) = getBoundaryMark(dofs_idx.at(i));
    //   }
    //   //i 为起点,j为终点,0-1:5,0-2:4,1-2:3。
    //   for(int i = 1; i >= 0; i--) {
    //     for(int j = 2; j > i; j--) {
    //       Point<2> new_point = (points.at(i) + points.at(j))/2.0;
    //       Dofs<2> new_dof(new_point,cglobal_index);
    //       dofs_idx.push_back(cglobal_index);
    //       new_dofsList.push_back(new_dof);
    //       int new_bndmark = (bndmarks.at(i) == 0 || bndmarks.at(j)==0) ? 0: bndmarks.at(i);
    //       if(new_bndmark != 0) {
    //         new_boundaryDofsList.push_back(new_dof);
    //       }
    //       // bndmarks.push_back(new_bndmark);
    //       pair<int,int> idx_bndmark = make_pair(cglobal_index,new_bndmark);
    //       new_indexonMesh.push_back(idx_bndmark);
    //       cglobal_index++;
    //     }
    //   }
    //   new_dofs_old_grid.push_back(dofs_idx);
    //   //构建新的边
    //
    // }
};

void EasyMesh::distributeDofs() {
    if(Order <=1 ) {
      return;
    }
    cout << "number of grids:" << nGrids << "\n";
    cout << "number of edges:" << nEdges << "\n";
    cout << "number of boundaryEdges:" << nBoundaryEdges << "\n";
    cout << "原自由度总数:" << nDofs << "\n";
    int npoints = nDofs;
    int new_nDofs = npoints + nEdges * (Order - 1) +
                    nGrids * ((Order - 1) * (Order - 2)) / 2;
    cout << "新自由度总数:" << new_nDofs << endl;
    nDofs = new_nDofs;
    dofsList.resize(nDofs);
    int numdofs_per_cell = ((Order + 1) * (Order + 2)) / 2;
    for (auto &indexs : indexonGrid) {
        indexs.resize(numdofs_per_cell);
    }
    indexonMesh.resize(nDofs);
    int global_index = npoints;
    for (int grid_idx = 0; grid_idx < gridsList.size(); ++grid_idx) {
        //cout << "deal with grid " << grid_idx << "\n";
        //取出自由度
        vector<int> dofs_idx = indexonGrid.at(grid_idx);
        vector<Dofs<2>> dofs;
        for (auto idx : dofs_idx) {
            dofs.push_back(dofsList.at(idx));
        }
        //一条线上自由度的边界标识由起点决定。
        vector<int> boundarymarks;
        vector<Point<2>> points;
        for (auto dof : dofs) {
            points.push_back(dof.getPoint());
            pair<int, int> index_bndmark = indexonMesh.at(dof.getGlobalIndex());
            int boundarymark = index_bndmark.second;
            boundarymarks.push_back(boundarymark);
        }
        //取出边，以在不同grid之间同步自由度
        vector<edgeIdx> edges_idx = edgeofGrid.at(grid_idx);
        assert(edges_idx.size() == 3);
        vector<vector<pair<elementIdx, edgeIdx>>> elements_idx;
        for (auto idx : edges_idx) {
            elements_idx.push_back(gridHasEdge.at(idx));
        }

        //分发边上的自由度
        //  2
        //  |\
        //  5—4
        //  | |\
        //  6—9—3
        //  | | |\
        //  0—7—8—1
        //  edge的编号为0,1,2
        //  edge0 : 1->2, edge1 : 2->0, edge2:*1->0*
        //  非常地奇怪
        //由于edge2的特殊性，先处理edge0,edge1;
        //cout << "deal with edge0 and edge1:"
        //     << "\n";
        int local_index = 3;
        for (int edge_idx = 0; edge_idx < 2; ++edge_idx) {
            int start_idx = edge_idx + 1;
            int end_idx = (start_idx == 2) ? 0 : (start_idx + 1);
            int start_mark = boundarymarks.at(start_idx);
            int end_mark = boundarymarks.at(end_idx);
            int boundarymark = (start_mark==0&&end_mark==0)?(0):(start_mark);
            Point<2> start = points.at(start_idx);
            Point<2> end = points.at(end_idx);
            for (int section = 1; section < Order; ++section) {
                int current_index = indexonGrid.at(grid_idx).at(local_index);
                if (current_index != 0) {
                    //已经在其他单元添加过，下一个自由度。
                    local_index++;
                    continue;
                }
                Point<2> new_point =
                    start + (end - start) * ((1.0 * section) / (1.0 * Order));
                Dofs<2> new_dof(new_point, global_index);
                dofsList.at(global_index) = new_dof;
                //cout << "edge0,1 local index is:" << local_index << "\n";
                // TODO:检查在local_index处是否以及有非零值，有的话不需要在添加，否则重复
                indexonGrid.at(grid_idx).at(local_index) = global_index;
                pair<int, int> index_bndmark =
                    make_pair(global_index, boundarymark);
                indexonMesh.at(global_index) = index_bndmark;
                if (boundarymark != 0) {
                    boundaryDofsList.push_back(new_dof);
                }
                //检查当前边是否属于另一个单元
                vector<pair<elementIdx, edgeIdx>> elements =
                    elements_idx.at(edge_idx);
                if (elements.size() > 1) {
                    for (auto ee : elements) {
                        if (ee.first != grid_idx) {
                            int element2_idx = ee.first;
                            int edge2_idx = ee.second;
                            int local2_index = indexBetweenElement(
                                edge_idx, local_index, edge2_idx);
                            indexonGrid.at(element2_idx).at(local2_index) =
                                global_index;
                        }
                    }
                }

                local_index++;
                global_index++;
            }
        }
        //处理edge2
        //cout << "deal with edge2:\n";
        int edge_idx = 2;
        int start_idx = 1;
        int end_idx = 0;
                  int start_mark = boundarymarks.at(start_idx);
            int end_mark = boundarymarks.at(end_idx);
            int boundarymark = (start_mark==0&&end_mark==0)?(0):(start_mark); 
        Point<2> start = points.at(start_idx);
        Point<2> end = points.at(end_idx);
        for (int section = Order - 1; section >= 1; --section) {
            int current_index = indexonGrid.at(grid_idx).at(local_index);
            if (current_index != 0) {
                //已经在其他单元添加过，下一个自由度。
                local_index++;
                continue;
            }
            Point<2> new_point =
                start + (end - start) * ((1.0 * section) / (1.0 * Order));
            Dofs<2> new_dof(new_point, global_index);

            dofsList.at(global_index) = new_dof;
            //cout << "edge2 local index is:" << local_index << "\n";
            indexonGrid.at(grid_idx).at(local_index) = global_index;
            pair<int, int> index_bndmark =
                make_pair(global_index, boundarymark);
            indexonMesh.at(global_index) = index_bndmark;
                if (boundarymark != 0) {
                    boundaryDofsList.push_back(new_dof);
                }
            //检查当前边是否属于另一个单元
            vector<pair<elementIdx, edgeIdx>> elements =
                elements_idx.at(edge_idx);
            if (elements.size() > 1) {
                for (auto ee : elements) {
                    if (ee.first != grid_idx) {
                        int element2_idx = ee.first;
                        int edge2_idx = ee.second;
                        int local2_index = indexBetweenElement(
                            edge_idx, local_index, edge2_idx);
                        indexonGrid.at(element2_idx).at(local2_index) =
                            global_index;
                    }
                }
            }

            local_index++;
            global_index++;
        }
        //分发内部自由度
        //cout << "deal with interior dofs:"
        //     << "\n";
        Point<2> d0 = points.at(0);
        Point<2> d01 = points.at(1) - d0;
        Point<2> d02 = points.at(2) - d0;
        for (int nume1 = ceil(Order / 2.0); nume1 >= 1; --nume1) {
            for (int nume2 = ceil(Order / 2.0); nume2 >= 1; --nume2) {
                if (nume1 + nume2 >= Order) {
                    continue;
                }
                Point<2> new_point = d0 +
                                     d01 * ((1.0 * nume1) / (1.0 * Order)) +
                                     d02 * ((1.0 * nume2) / (1.0 * Order));
                Dofs<2> new_dof(new_point, global_index);

                dofsList.at(global_index) = new_dof;
                //cout << "interior local index is:" << local_index << "\n";
                indexonGrid.at(grid_idx).at(local_index) = global_index;
                pair<int, int> index_bndmark =
                    make_pair(global_index, 0);
                indexonMesh.at(global_index) = index_bndmark;
                local_index++;
                global_index++;
            }
        }
    }
};

void EasyMesh::inputDofData() {
    ifstream fin;
    fin.open(meshName + ".n");
    string line;
    stringstream ss;
    getline(fin, line);
    ss.clear();
    ss.str(line);
    ss >> nDofs >> nGrids >> nEdges;
    gridHasEdge.resize(nEdges);
    int index, bndflag;
    double x, y;
    for (int i = 0; i < nDofs; i++) {
        getline(fin, line);
        ss.clear();
        ss.str(line);
        ss >> index >> x >> y >> bndflag;

        Point<2> temp_point({x, y});
        Dofs<2> temp_dof(temp_point, index);
        indexonMesh.push_back({index, bndflag});
        dofsList.push_back(temp_dof);
        if (bndflag != 0)
            boundaryDofsList.push_back(temp_dof);
    }
}

void EasyMesh::inputEdgeData() {
    ifstream fin;
    fin.open(meshName + ".s");
    string line;
    stringstream ss;
    getline(fin, line);
    ss.clear();
    ss.str(line);
    int a, d0, d1, e0, e1, bndflag;
    for (int i = 0; i < nEdges; i++) {
        getline(fin, line);
        ss.clear();
        ss.str(line);
        ss >> a >> d0 >> d1 >> e0 >> e1 >> bndflag;
        vector<int> temp_edge(2);
        temp_edge[0] = d0;
        temp_edge[1] = d1;
        edgesList.push_back(temp_edge);
        if (bndflag != 0) {
            boundaryEdgesList.push_back(temp_edge);
        }
    }
    nBoundaryEdges = boundaryEdgesList.size();
}

void EasyMesh::inputGridData() {
    ifstream fin;
    fin.open(meshName + ".e");
    string line;
    stringstream ss;
    getline(fin, line);

    int a, d0, d1, d2, e0, e1, e2, s0, s1, s2;
    for (int i = 0; i < nGrids; i++) {
        getline(fin, line);
        ss.clear();
        ss.str(line);
        ss >> a >> d0 >> d1 >> d2 >> e0 >> e1 >> e2 >> s0 >> s1 >> s2;
        vector<int> indexs = {d0, d1, d2};
        indexonGrid.push_back(indexs);
        vector<Point<2>> points = {dofsList[d0].getPoint(),
                                   dofsList[d1].getPoint(),
                                   dofsList[d2].getPoint()};
        Grid<2> grid(points);
        gridsList.push_back(grid);
        //每个grid的edge;
        vector<edgeIdx> edges = {s0, s1, s2};
        edgeofGrid.push_back(edges);
        //每个edge对应的grid;
        int eidx = i;
        auto p0 = make_pair(eidx, 0);
        gridHasEdge.at(s0).push_back(p0);
        auto p1 = make_pair(eidx, 1);
        gridHasEdge.at(s1).push_back(p1);
        auto p2 = make_pair(eidx, 2);
        gridHasEdge.at(s2).push_back(p2);
    }
};

Grid<2> &EasyMesh::getGrid(int index) { return gridsList.at(index); };

vector<Grid<2>> &EasyMesh::getGrid() { return gridsList; };

vector<Grid<2>> EasyMesh::getBoundaryGrid() {
    // why we need to getBoundaryGrid? we already have getBoundaryDofs
    // TODO:this function is not used,we can delete it.
    return gridsList;
};

set<Dofs<2>> EasyMesh::getBoundaryDofs() {
    // use set to avoid duplicate dofs
    // unnecessary in unstructed mesh
    set<Dofs<2>> dofsSet(boundaryDofsList.begin(), boundaryDofsList.end());
    return dofsSet;
};

vector<pair<int,int> > EasyMesh::getIndexofGrid(int i)
{
  vector<pair<int,int> > res;
  vector<int> dofsindex = indexonGrid.at(i);
  for (int& index : dofsindex)
  {
    auto boundarymark = find_if(indexonMesh.begin(), indexonMesh.end(), [index](const pair<int, int>& p) {
        return p.first == index;
    });
    if (boundarymark != indexonMesh.end()) {
      res.emplace_back(index,(*boundarymark).second);
    } else {
      std::cerr << "no such index in indexonMesh" << std::endl;
    }
  }
  return res;
};
#else
#endif
